A Trajectory Analysis using the ICESat2 Satellite

Overview of the steps involved:

  1. Get very high-fidelity orbit trajectory information of the primary satellite:

    • This can be done using a POD-reduced dynamics run in which all empirical accelerations and other parameters are used to compensate for any mismodelled forces.

    • This has already been done for our purposes.

  2. Do a run of GEODYN in which the orbit trajectory IS the tracking datatype (PCE).

    • Use this run type to do all density model assessments

    • The residuals will be the difference from the truth.

  3. Methods for assessing the density models (using trajectory run type):

  • General fit (residuals and RMS)

  • Arc overlap

    • If any arcs overlap, look at how well different density models provide consistency between overlapping arcs

  • Test of prediction - Do a fit from t1 to t2 and then predict to some t3. If the density model is better, the difference between predicted orbit (t2 to t3) and the precise trajectory (PCE data) will contain how well the model is doing.

Read GEODYN Output using PygeodynRead functionality

Get MSIS2 Data

[1]:

# Clear the memory of the notebook because the saved variables get so large:
# %reset -f out


[2]:

%load_ext autoreload
%autoreload 2

import sys
sys.path.insert(0, '/data/geodyn_proj/pygeodyn/pygeodyn_develop/')
from PYGEODYN import Pygeodyn

#------ A dictionary containing the run parameters ------
run_params1 = {}
run_params1['arc']               =  ['2018.313',
                                     '2018.314',
                                     '2018.315',
#                                      '2018.316',
#                                      '2018.317',
                                    ]
run_params1['satellite']         =  'icesat2'
run_params1['den_model']         =  'msis2'
run_params1['SpecialRun_name']   =  '_TrajAnalysis'
run_params1['verbose']           =  False
run_params1['action']            =  'read'
run_params1['request_data'] = ['AdjustedParams',
            #                    'Trajectory_xyz',
                               'Trajectory_orbfil',
                               'Density',
                               'Residuals_obs',
                               'Residuals_summary',
                               'Statistics',
                                ]

Obj_Geodyn1 = Pygeodyn(run_params1)
Obj_Geodyn1.getData()


                      ......... READING GEODYN output
     Loading ... icesat2_2018313_54hr.msis2
     Loading ... icesat2_2018314_54hr.msis2
     Loading ... icesat2_2018315_54hr.msis2

Get MSIS00 Data

[3]:

%load_ext autoreload
%autoreload 2

from PYGEODYN import Pygeodyn
import copy


run_params2 = copy.deepcopy(run_params1)

run_params2['den_model']         =  'msis00'

Obj_Geodyn2 = Pygeodyn(run_params2)
Obj_Geodyn2.getData()


The autoreload extension is already loaded. To reload it, use:
  %reload_ext autoreload
                      ......... READING GEODYN output
     Loading ... icesat2_2018313_54hr.msis00
     Loading ... icesat2_2018314_54hr.msis00
     Loading ... icesat2_2018315_54hr.msis00

Get MSIS86 Data

[4]:

%load_ext autoreload
%autoreload 2

from PYGEODYN import Pygeodyn
# import copy


run_params3 = copy.deepcopy(run_params1)
run_params3['den_model']         =  'msis86'

Obj_Geodyn3 = Pygeodyn(run_params3)
Obj_Geodyn3.getData()


The autoreload extension is already loaded. To reload it, use:
  %reload_ext autoreload
                      ......... READING GEODYN output
     Loading ... icesat2_2018313_54hr.msis86
     Loading ... icesat2_2018314_54hr.msis86
     Loading ... icesat2_2018315_54hr.msis86

Analyze output

Our analysis has the following output products:

  1. Residuals of the POD across many arcs

  2. RMS of fit of the POD across many arcs

  3. Adjustment of the drag coefficient (drag acceleration to compensate for inaccuracies in the density model.)

  4. Check of consistency (in the residuals) across overlapping arc times

  5. RMS of the overlapping residual difference (this removes the PCE contribution) Using the Orbfil:

  6. From ORBFIL grab the overlapping ephemeris and difference the two models. Compare this against the same of prediction

    • How well does the predicted time period match up with the determined ephemeris (see this in the resids of the two).

  7. Calculate and plot the radial component of the trajectory

[5]:

import plotly.graph_objects as go
from plotly.offline import plot, iplot
from plotly.subplots import make_subplots
import plotly.express as px


config = dict({
                'displayModeBar': True,
                'responsive': False,
#                 'staticPlot': True,
                'displaylogo': False,
                'showTips': False,
                })

1. Residuals of the POD across many arcs

Resids = PCE - POD Trajectory

[6]:

%load_ext autoreload
%autoreload 2

from PYGEODYNAnalysis_icesat2PCEtrajectory import plot_residual_meas_summary
from PYGEODYNAnalysis_icesat2PCEtrajectory import rms_summary_table


Obj_list = [Obj_Geodyn1,Obj_Geodyn2,Obj_Geodyn3,]
rms_summary_table(Obj_list)



fig = make_subplots(rows=2, cols=1,
     subplot_titles=(["Mean Residuals per Arc", 'RMS of Fit per Arc']),
     vertical_spacing = 0.1)
fig = plot_residual_meas_summary(fig, Obj_Geodyn2, 0)
fig = plot_residual_meas_summary(fig, Obj_Geodyn3, 1)
fig = plot_residual_meas_summary(fig, Obj_Geodyn1, 2)
fig.show(config=config)



The autoreload extension is already loaded. To reload it, use:
  %reload_ext autoreload
+———————————————————+—————————————————————————+————————————————+
|                     Summary Across all Arcs                  |
+———————————————————+—————————————————————————+————————————————+
+   Density Model   +   Mean Residual (cm)    +   RMS of Fit   +
+-------------------+-------------------------+----------------+
+       msis2       +       3.13333e+01       +  1.88767e-01   +
+       msis00      +       5.00000e+01       +  1.95100e-01   +
+       msis86      +       2.63000e+02       +  3.73433e-01   +
+———————————————————+—————————————————————————+————————————————+

[7]:

%load_ext autoreload
%autoreload 2

from PYGEODYNAnalysis_icesat2PCEtrajectory import plot_residuals_observed



fig = make_subplots(rows=3, cols=1,
            subplot_titles=(['X', 'Y', 'Z']),
            vertical_spacing = 0.1,
                       )

fig = plot_residuals_observed(fig, Obj_Geodyn2, 0)
fig = plot_residuals_observed(fig, Obj_Geodyn3, 1)
fig = plot_residuals_observed(fig, Obj_Geodyn1, 2)

fig.update_layout(title="Observation Residuals (PCE - Observed , T.O.R.)")

fig.show(config=config)


The autoreload extension is already loaded. To reload it, use:
  %reload_ext autoreload

[8]:

%load_ext autoreload
%autoreload 2
from PYGEODYNAnalysis_icesat2PCEtrajectory import plot_cd_and_percdiff_from_apriori


fig = make_subplots(
    rows=2, cols=1,
    subplot_titles=(["Timeseries of Cd", "Percent difference from a priori (Cd=2.2)"]),
    vertical_spacing = 0.08,
    )
fig = plot_cd_and_percdiff_from_apriori(fig,  Obj_Geodyn2, 0)
fig = plot_cd_and_percdiff_from_apriori(fig,  Obj_Geodyn3, 1)
fig = plot_cd_and_percdiff_from_apriori(fig,  Obj_Geodyn1, 2)


fig.show(config=config)

The autoreload extension is already loaded. To reload it, use:
  %reload_ext autoreload

[9]:

%load_ext autoreload
%autoreload 2

from PYGEODYNAnalysis_icesat2PCEtrajectory import plot_ScaleDensity_with_CdScaleFactor__2

fig = make_subplots(rows=2, cols=1,
                subplot_titles=(["Model Ouptut Density", "Model Density * Cd Scaling Factor"]),
                shared_yaxes=True,
                vertical_spacing = 0.1,
                specs=[
                [{"secondary_y": False}],
                [{"secondary_y": False}], ])


fig = plot_ScaleDensity_with_CdScaleFactor__2(fig,  Obj_Geodyn2, 0, 200)
fig = plot_ScaleDensity_with_CdScaleFactor__2(fig,  Obj_Geodyn3, 1, 200)
fig = plot_ScaleDensity_with_CdScaleFactor__2(fig,  Obj_Geodyn1, 2, 200)

# min_y = 1*1e-16
# max_y = 9*1e-12
# fig.update_yaxes(range=[min_y, max_y], row=1, col=1)
# fig.update_yaxes(range=[min_y, max_y], row=2, col=1)

fig.show(config=config)


The autoreload extension is already loaded. To reload it, use:
  %reload_ext autoreload

Looking at the arc overlap time:

We want to show the residuals in the overlap time with the PCE data subtracted away.

[10]:

####  ARC_OVERLAP_ObsResids_XYZ

%load_ext autoreload
%autoreload 2

from PYGEODYNAnalysis_icesat2PCEtrajectory import ARCOVERLAP_2arcs_ObsResids_XYZ

fig = make_subplots(rows=3, cols=1,
            subplot_titles=(['X', 'Y', 'Z']),
            vertical_spacing = 0.1,
            specs=[ [{"secondary_y": True }],
                    [{"secondary_y": True }],
                    [{"secondary_y": True }], ],)

arc1 = '2018.314'  # '2018.314'
arc2 = '2018.315'

fig = ARCOVERLAP_2arcs_ObsResids_XYZ(fig, Obj_Geodyn2, 0, arc1, arc2)
fig = ARCOVERLAP_2arcs_ObsResids_XYZ(fig, Obj_Geodyn3, 1, arc1, arc2)
fig = ARCOVERLAP_2arcs_ObsResids_XYZ(fig, Obj_Geodyn1, 2, arc1, arc2)

fig.show(config=config)


The autoreload extension is already loaded. To reload it, use:
  %reload_ext autoreload

PCE Data and the Orbit File:

Residual Component Trajectory:

Convert the Interial XYZ coordinates to the satellite coordinate system (RSW), then plot the radial component.

Starting Systems: - PCE data - J2000 Coordinate System - Inertial satellite State Vector: \([x, y, z, \dot{x}, \dot{y}, \dot{z}]\) (m) - ORBFIL data - Mean of year Coordinate System - Inertial satellite State Vector: \([x, y, z, \dot{x}, \dot{y}, \dot{z}]\) (m)

Convert from ``XYZ`` to ``RSW``

From Vallado pg. 164:

[11]:


%load_ext autoreload
%autoreload 2

from PYGEODYNAnalysis_icesat2PCEtrajectory import ARCOVERLAP_2arcs_ObsResids_RSW_radial

fig = make_subplots(rows=2, cols=1,
            subplot_titles=(['Radial Component', 'Residual (PCE-ORBFIL)']),
            vertical_spacing = 0.2,
            specs=[ [{"secondary_y": False }],
                    [{"secondary_y": False }]],)

arc1 = '2018.314'
arc2 = '2018.315'

fig = ARCOVERLAP_2arcs_ObsResids_RSW_radial(fig, Obj_Geodyn2, 0, arc1, arc2)
fig = ARCOVERLAP_2arcs_ObsResids_RSW_radial(fig, Obj_Geodyn3, 1, arc1, arc2)
fig = ARCOVERLAP_2arcs_ObsResids_RSW_radial(fig, Obj_Geodyn1, 2, arc1, arc2)

fig.show(config=config)


The autoreload extension is already loaded. To reload it, use:
  %reload_ext autoreload

[ ]:


%load_ext autoreload
%autoreload 2

from PYGEODYNAnalysis_icesat2PCEtrajectory import ARCOVERLAP_2arcs_ObsResids_NTW_intrack

fig = make_subplots(rows=2, cols=1,
            subplot_titles=(['In-Track Component', 'Residual (PCE-ORBFIL)']),
            vertical_spacing = 0.2,
            specs=[ [{"secondary_y": False }],
                    [{"secondary_y": False }]],)

arc1 = '2018.314'
arc2 = '2018.314'

fig = ARCOVERLAP_2arcs_ObsResids_NTW_intrack(fig, Obj_Geodyn2, 0, arc1, arc2)
fig = ARCOVERLAP_2arcs_ObsResids_NTW_intrack(fig, Obj_Geodyn3, 1, arc1, arc2)
fig = ARCOVERLAP_2arcs_ObsResids_NTW_intrack(fig, Obj_Geodyn1, 2, arc1, arc2)

fig.show(config=config)


The autoreload extension is already loaded. To reload it, use:
  %reload_ext autoreload

[ ]: